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We use functional Bethe Ansatz equations to calculate the cumulants of the total current in the 
partially asymmetric exclusion process. We recover known formulas for the first two cumulants 
(mean value of the current and diffusion constant) and obtain an explicit finite size formula for 
the third cumulant. The expression for the third cumulant takes a simple integral form in the 
limit where the asymmetry scales as the inverse of the square root of the size of the system, which 
corresponds to a natural separation between weak and strong asymmetry. 

PACS numbers: 05-40.-a; 05-60.-k 

Keywords: ASEP, functional Bethe Ansatz, large deviations 



I. INTRODUCTION 



The asymmetric simple exclusion process (ASEP) is one of the most simple examples of a stochastic interacting 
particles model with an out of equilibrium stationary state. It features classical particles hoping on a lattice and 
interacting through hard core exclusion. It can be seen as a growth model with deposition and evaporation of 
particles and is thus a discrete version of a system evolving by the Kardar-Parisi-Zhang (KPZ) equation. The one 
dimensional model has been much studied in the past [TJ [21 |3j HI |5] . It is known to be exactly solvable through various 
methods including the Bethe Ansatz [B] [7j , which uses the underlying integrability of the Markov matrix of the model 
(similar to the Hamiltonian of an XXZ spin chain with twisted boundary conditions) , and the matrix product Ansatz 
method [SjH], which consists in writing the stationary state eigenvector as a trace of product of matrices. 

The fluctuations of the current have received much attention [lOl [TTJ [12] . In the long time limit, the system reaches 
a stationary state, independent of the initial configuration for finite size systems. A manifestation of the fact the 
system is not at equilibrium in the stationary state is the existence of a current of particles flowing through the 
system. In this stationary state, the mean value J of the current is simply proportional to the asymmetry of the rates 
at which the particles hop to the right and to the left, which can be thought as a driving field. The fiuctuations of 
the current describe how the current differs from its mean value, that is the probability of finding a current j different 
from J. But finding the whole probability density function can be difficult. Instead, one can study its first cumulants. 
The second cumulant, related to the diffusion constant, describes the width of the probability density whereas the 
third cumulant represents the asymmetry and the non-gaussianity of the probability density. The fourth cumulant 
describes how much the peak of the distribution is sharp. 

The ffuctuations of the current have been studied for different variants of the ASEP: finite size open lattice [T^], 
systems with second class particles [li|[l3], infinite lattice [14l[15],... For the totally asymmetric exclusion process, in 
which the particles only hop in one direction (TASEP), on a finite size lattice with periodic boundary conditions, all the 
cumulants of the probability distribution of the current have been calculated by Bethe Ansatz, using a simplification 
of the Bethe Ansatz equations that occur only in this case. For the partially asymmetric exclusion process on a ring, 
a particular scaling limit of the probability distribution function of the current was calculated using a thermodynamic 
limit of the Bethe Ansatz equations |16j . However, only the first two cumulants were known for finite size systems 
pn [T5] . In the present paper we introduce a method that allows us to calculate exact expressions for the first 
cumulants. We use it to get the diffusion constant with very few calculations, and obtain an explicit formula (58 1 for 
the third cumulant. We also derive a rather simple integral formula (69 1 for the third cumulant in the scaling limit 
where the asymmetry goes to zero as the inverse of the square root of the size of the system. This formula is similar 
to the one obtained in [17 for the diffusion constant in the same scaling limit. 

The structure of the paper is as follows: In section [ll] we recall the functional Bethe Ansatz equations used in 
[TSl to calculate the two first cumulants of the total current. In section [ill we show how to solve perturbatively 



these equations, giving a way to calculate explicitly the cumulants one after the other. We obtain easily the two 
first cumulants (mean value of the current and diffusion constant). For the diffusion constant, our derivation is much 
simpler than in |17j and |18j . In section |V] after a more involved calculation, we obtain an explicit expression for 
the third cumulant. We show that for the cases of either the symmetric exclusion process (SSEP, the particle hop on 
both sides with the same rate) or TASEP it reduces to known formulas. We finally write the third cumulant formula 
in an integral form in the scaling limit where the asymmetry goes to zero as the inverse of the square root of the size 
of the system. 
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II. BETHE ANSATZ FOR THE FLUCTUATIONS OF THE CURRENT 



We consider in this paper the partially asymmetric simple exclusion process (PASEP) on a one dimensional lattice 
with periodic boundary conditions. Each of the L sites can be empty or occupied by at most one of the n particles. The 

configuration space of the model has then dimension = ^ ^ ^ • this configuration space, we define the following 

Markov process: each particle can hop to the right with rate p (that is, with probability pdt for an infinitesimal time 
interval dt) and to the left with rate q provided that the target site is empty. Otherwise, the particle cannot move. 
This dynamics can be encoded in the Markov matrix M = Mq + pMi + qM^i, whose diagonal part Mq contains the 
exit rates from each configuration, and Mi (resp. M_i) the incoming rates obtained by moving one particle to the 
right (resp. to the left). 

Let Yt be the total distance covered by all the particles between time and t. Yf is thus the integrated current 
between time and t. It has been shown [TSl |TS] that the generating function of the cumulants of Yt in the 
long time limit can be obtained from the diagonalization of a deformation of the Markov matrix of the system: 
M(7) = Mq + pe^Mi + qe^'^M^i, introducing the parameter 7 which can be seen as a fugacity associated with the 
leaps of the particles. More precisely, the eigenvalue £'(7) of M{j) corresponding to the stationary state of the system, 
that is its eigenvalue with maximal real part, is given by 

Ei,)^l^J^^^J,+ ^,^+f,^ + -- (1) 
J is the total current and A is the diffusion constant. 

J = hm ^ (2) 



t^oo t 

J?) ~ {Yt 

t — >oo 



lim (3) 



E, = lim 3(r.)(r.-) + 2(yj3 

t^oo t 

The diagonalization of -M(7) can be performed using the Bethe Ansatz, because of the underlying integrability of the 
model. We are now going to recall the functional equations derived in [18 for the calculation of E^-f). Defining the 
asymmetry parameter x by 

(5) 

P 

the Bethe equations of the system are 



e 



The Bethe roots yi I < i < n depend on both x and 7. They can be used to write the corresponding eigenvector as a 
sum over the n\ permutations of the particles. Introducing the polynomial 

n 

Q{t) = l[{t-y^) (7) 

1=1 

whose zeros are the Bethe roots, it can be shown [18] that the Bethe equations (|6| can be rewritten as functional 
equations 

Qit)R{t) = e^^(l - t)^Q{xt) + (1 - xt)^x"'Q{t/x), (8) 

where Q is a polynomial of degree n (with a coefficient of highest degree equal to 1) and i? is a polynomial of degree L. 
Both Q and R must be determined by solving equation (|8|. This equation is known as Baxter's scalar TQ equation 
[5ni[5TJ[51[23]. Equation ^ sets the value of R{0) 

R(0) = + e^T (9) 

and the behavior of R(t) when t — > 00 

R{t) - (x^ + x''e^'^){-l)^t^ when t -> 00. (10) 
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When 7 = 0, the solution of the Bethe equations ([6| corresponding to the stationary state is ?/i = for all i. 
Equivalently, the corresponding solution of dsl) is Q{t) = t" and R{t) = a;"(l — t)^ + (1 — xt)^. Expanding Q{t) and 
R{t) near 7 = 0, we get 



Q(t) = r + o (7) 



and 



R{t) = - i)^ + (1 - .Tf)^ + O (7) ■ 
The stationary state has another interesting property: it is a state with zero momentum, which leads to 

x"Q{l/x) = e"^g(l) 

or, in terms of the polynomial R 



(11) 
(12) 

(13) 
(14) 



The latter equation does not provide additional information beside (12 1, but it is useful to make some calculations 
easier. 

Finally, with Q and R solution of (l8|, the eigenvalue of Af(7) corresponding to the stationary state is given by [TS] 



^M^n- ^/^Q'(l) i Q'(i/^) A 

p ^ ^{q{1) xQ{1/x)) 



= —Lx — {1 — x) 



R{1)- 



(15) 



Equation ([8]) could be solved perturbatively near 7 = right now. Inserting the solution into (15 1, it would give us 
the first cumulants of the current. This was done in [18j up to order 2. Here, we will rewrite the previous equations 
in a different way before making the perturbative expansion near 7 = 0. It will make the calculations much simpler 
than in [18], and allow us to calculate the third cumulant. We divide both sides of the functional equation ^ by 
{l-t)^{l-xt)^Q{t) 



x"-e^'^ 



m ^ A{t) 1 

(1 - t)L{l - xt)L (1 - t)L (1 - xt)L A{xt) ' 



where we defined 



A{t) = x' 



Q{t) ■ 



(16) 



(17) 



From this definition of A{t), and using (13) and the fact that Q is of degree n, we know the value of A{t) in t = 0, 
t = 1 and in the limit i — > 00 



A{0) = x", A{1) = e"^ 



and 



hm A{t) = 1 

t — *oo 



From (11 1, we also know that 



A(t) = 1 + 0(7). 



In terms of A{t), the expression (151 for the eigenvalue rewrites, using the value of ^(1) (18) 

= -(1 -a;)e-"^yl'(l). 



(18) 



(19) 



(20) 



III. PERTURBATIVE SOLUTION OF THE FUNCTIONAL EQUATION 



In this section, we will first show how to eliminate R{t) from equation (16), leaving us with a closed equation for 
A(t). We will see how we can solve this equation perturbatively near 7 = 0. Then, we will reformulate this solution 
to make the explicit calculation of the cumulants easier. 



4 



A. Perturbative solution for A{t) 



From (19), we can write the expansion of A{t) near 7 = as 

00 

A(t) = l + ^^,.(t)7* 



(21) 



But Q{t) is normalized such that Q{t) — is a polynomial of degree n—1. Thus, from ( 17 1, the Ak{t) are polynomials 
in 1/i of degree kn without constant term {Ak{t) when t 00). 



From now on, we will look at a perturbative solution of (16 1 for small 7. R{t), which is a polynomial in t, can be 



seen as a formal series in t and 7 with only nonnegative powers in t. The l.h.s. of equation ( 16 ) is then also a formal 
series in t and 7 with only nonnegative powers in t. But we also know that, at each (nonzero) order in 7, A{t) has 
only negative powers in t so that the r.h.s. of equation (16 1 has both negative and nonnegative powers in t. Thus, 



equation ( 16 ) simply means that the negative powers in t in equation ( 16 1 cancel. We will write this as 



1 



(1 -xt)L A{xt) 



(22) 



which will mean: at each order in 7, the l.h.s. of ( |22| is finite when t 
(16 1. We will see that equation (22 1 still determines A{t) uniquely. 



0. We have eliminated R(t) from equation 



From now on, every expansion in powers of t will have to be understood as an expansion of a formal series in 
powers of 7 followed, at each order in 7, by an expansion in powers of t. In the following, we will use the notation 
[/(^)](/c) = [/](fe) to refer to the coefficient of the term t*^ in the expansion of the formal series f{t). We also introduce 
the notation [/(i:)](_) which will represent the negative powers of / in t (that is terms t'^ with fc < 0). On the contrary, 
[f{t)]i+) will represent the nonnegative powers of / in t (terms t'^ with fc > 0). 

To solve equation ( 22 1 , we will write it in the slightly more complicated form 



{1-ty 



(1 - xty 



A{xt) 



A{xt) 



Ax being the operator which acts on an arbitrary function u{t) = ^ [^Jfj-^i*^ as 

fcez 



{Axu){t) = u{t) - x'^u{xt) 



E 



(1 



(23) 



(24) 



As Ax gives when applied to t the equation 



{Axu){t) = v{t) 



can be formally solved as 



k^—n 



v]{k)t'' _ b 

j,n-\-k ^r, 



(25) 



(26) 



for any formal series v which has no term t 6 is an arbitrary coefficient which is not constrained by equation (25 1 



We can use this to invert A^, in (23 1. Recalling that A{t) 1 when t ^ 00 and that A(t) — 1 has only negative 
powers in t, we obtain 



^(i) = l-[(l-<)^g(i)](_)-6 



with 



m 



A: 



A{xt) 



(-) 



(27) 



(28) 



The term containing b is not constrained by equation ( 23 1 . The condition necessary to invert A^, in equation ( 23 1 



which is that the r.h.s. of (23 1 must not contain a term t implies that g{t) does not contain a term t 
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Written like that, equation (27 1 contains A on both sides (through g on the r.h.s.). However, equation (27 1 gives a 
recursive solution for A{t) order by order in powers of 7. Indeed, if we take the term of order 7*^ in this equation, the 
l.h.s. depends only on Ai^{t) whereas the r.h.s. depends only on the Aj{t) with j < k: the Ai^{t) cancels out at order 
k in the expansion in powers of 7 of 

A{xt) 



A{xt) 



fc=i 1 



1 + {e^^ - 1) 



(29) 



^ Ak{xt)-f'' 

k=l 



Thus, equation (27) allows to recursively compute all the orders in 7 starting from (19 1. The parameter b in (27 1 can 
be set from the value (flSl) of A{1). 



Using equations (271), il9h and (18), we can do the calculation for the first cumulants. But wc will first turn these 



equations into another form which will make the calculations a little simpler. 

We must emphasize that our method involves only algebraic manipulations of formal series in two parameters, 
namely 7 and t. The Bethe roots 2/^(7) (which are the zeros of the polynomial Q{t) and the poles of the rational 
function A{t)) have completely disappeared since we did the 7 expansion first: the Ak{t) have only a pole of order kn 
in t = 0. We do not have to follow the 2/^(7) as a function of 7: we only use some algebraic properties of A{t). 



B. A simpler formulation for the perturbative solution 



We will now eliminate A(t) from our recursive equations, and wor k on the auxiliary quantity g{t). Noticing that 
A~^{xt) — 1 has also only negative powers in t, wc find from equations (22) and (27) that x"e^'^A'~^{xt) can be written 
in terms of g(t) as 

[{l-xt)^g{t)]i-) + b^^^''''^'^ 



(-) 



(30) 



(31) 



We can absorb the b terms of A{t) and x"e^'^A~^{xt) into g, defining 

g{t)^g{t) + l^- 

b is thus the coefficient of i^" in g{t) as g{t) has no term In the following, we will both need using g{t) and g{t). 
The expressions for A{t) and x'^e'"'^ A~^{xt) become 

A(t) = !-[(!- i)^5(t)](-) (32) 
j^^e^-rA-^xt) = x''e^^ + [(1 - xt)^g{t)]^_y (33) 
From (28), at order fc in 7 g{t) is a polynomial in 1/t of degree kn without constant term, as are the Ak{t). We also 

(34) 



see that g{t) is equal to when 7 = (as A{t) = 1 when 7 = 0), giving the expansion 



k = l 



To find a closed equation for g{t), we eliminate A{t) between (32) and (33), using {A{t)){x"e^'^ A ^(t)) = x'^e^'^: 

[(1 - t)^g{t/x)\.) - x-e^'<[{l ~ t)^g{t)](.) - [(1 - t)^g(i)](-)[(l - t)^5(t/x)](-) = 0, (35) 
which can be rewritten as 

^^"^^^ ](-)[(l- 

As g(t) is of order 7, this allows us to solve order by order in powers of 7. Wc can simplify the r.h.s. of the previous 
equation, writing 



,itM - .^e^Mt) ^ [(l-^)-.(^)](-)Kl -t)^.(V.)](-) ^ ^ (,o) 



(36) 



[(l-t)^g(t)](-)[(l-^)^g(t/x)](_) 
(l-i)i 

^ ((1 - t)^g{t) [(1 - t)^.g(t)](+)) ((1 - t)^git/x) - [(1 - t)^g{t/x)]^+^) 

= (1 - t)^g{t)git/x) - 5(t/x)[(l - t)^gi%+) - g(i)[(l - t)^git/x)](^+^ + O 

= -(1 - t)^g{t)git/x) + g{t/x)[{l - t)^<?(i)](_) + 5(i)[(l - t)^ g{t/x)]^_^ + O {t°) 



(37) 
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and we obtain 

g{tlx) - x-e^-'g{t) = -(1 - t)^g{t)g{t/x) + - t)^ g{t/ x)]^^) + g{t/x)[{l - + O (38) 

or, at order r in 7, multiplying thyx 



(Aa;5r)(t) = a;" ^ -^gr-k{xt) 



(39) 



r-l 



^ ((1 - xt)'^gr-k{xt)gk{t) - gk{tm - xi) V-fc(2;t)](-) - gr^k{xt)[{l - x<) + O (f) , 



k=l 



where we used once again the operator (24 1. Using the formal inversion formula (26) for A^^, we see that the 
previous equation for g{t) can be solved if 



r j^k 



E 

fc=i 



x'^—gr-^kixt) - {1 ~ xt) gr-k{xt)gk{t) + gk{t)[{l ~ xt) gr-k{xt)](^) + gr-k{xt)[{l ~ xt) gk{t)](-) 



= 0, 



(-") 



(40) 



which is a consequence of the fact that equation (39) holds: it must be true if ([s]) has a solution verifying ( [Tl] ) and 
(fT2l. We obtain 



r-l (r-k)n 







k=l 1=1 



-I r^fc 
Id 



gr~k{xt) 



(-1) 



(41) 



EE 

k=l 1 = 1 



[(1 - xt)^gr-k{xt)gk{t) - .gfc(t)[(l - a;t)^g^_fc(a;t)](_) - gr-k{xt)[{l - xt)^ gk{t)](^)] 

(1 - .T"-')*' 



and 



(42) 



The coefficient 5r is the term of order 7'' in b. Again, the value (18) of A{1) sets br- We note that this formula does 



not require to divide by a formal series in 7, unlike (27) where we had to divide by A(xt). The nonlinearity of the 



recurrence formula is then simpler here; it consists only in a product of two series in 7. 



IV. EXPLICIT CALCULATIONS FOR THE TWO FIRST CUMULANTS OF THE CURRENT 

In this section, after expressing the generating function of the cumulants of the current in terms of the perturbative 
solution described in the previous section, we will calculate explicitly the two first cumulants. 



A. Expression of the stationary state eigenvalue 

We will now express the stationary state eigenvalue using the function g{t). As g{t) has only negative powers in t, 
expression (32) for A{t) in terms oi g{t) can be written as 



i=l ]=0 



^(^)-i-EE[5](-)( ^ 



J 



(43) 



Using the binomial coefficients formulas (Al I and (A2|, we get for the eigenvalue 

L(i-l)^-'^^(-^')- 



£^(7) 1-x 



j=i 



(44) 
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This expression for £^(7) is of order 7, as is g{t), which means E{'y = 0) = 0, as expected for the stationary state. 

The last equation allows us to calculate the eigenvalue at any order in 7 if we know g{t) at the corresponding order. 
g{t) is obtained by the recurrence equation (41) in which br- is set using the known value (18 1 of ^(1). Thus, in the 
calculation of the eigenvalue br can be eliminated from (44 1 where it appears through g{ty^ g{t) + If we want 



to calculate the eigenvalue at order 3, we will not need 63 (but we will still need 61 and 62 as the recurrence equation 
(41 1 involves g{t) and not only g{t)). Using the value (18 1 of ^(1), 



which gives for b = ^rl^ 



with 



(45) 



b e"'' - 1 ^ f L 



y n f \ J I n L ^ ■" 



(46) 



^^(-1)" 



We finally get for the eigenvalue, putting (44 1 and (46 1 together 



(47) 



(48) 



We will now use (41 1, (46 1 and (48 1 to calculate £'(7) up to order 2 in 7, recovering known results for the mean value 



of the current and for the diffusion constant. 



B. Calculation of the mean value of the current 

From ([T|), the mean value J of the current is given by 



(49) 



We only need E(j) at the first order in 7. In this case, the two complicated terms with the sums in (41 1 do not 
contribute and gi{t) = 0. Using (48 1, the eigenvalue reads 



p L — 1 



We did not need 61 to calculate J, but we will need it for the next orders. Using (46 1 we find bi — v and 

V 



9i{t) 



(50) 



(51) 



Calculation of the diffusion constant 



From ([T]), the difi^usion constant A is given by 



A 

P 



(52) 



|7=0 
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The recurrence ( 41 1 gives 



92 



L 

2n-k 



i-ty 



fe=l 



k — 7 



2n 
fc=n+l 



1 - X 



k — r 



(53) 



The eigenvalue is, using (48) and (A4| 



^ fe>o 



71 + fc 



L 

n - k I 1 



1 — a;'' 



(54) 



The derivation of this expression is much simpler than the previous derivation using the Bethe Ansatz [18'. It uses 
a little more formalism, but nearly no calculation is needed, contrary to the previous derivation in which many 
unexpected simplifications occurred in the end of the calculation. 



The condition ( 40 1 necessary for the consistency of equation ( 39 1 is easily checked 

ia;"[5i(i)](-n) - [(1 - i)^5i(t)3i(Va^)](-n) = 0. 



(55) 



As for the first order, we did not need 62 to calculate the diffusion constant, but we will need it for the third order. 



We set 62 using equation (46 1 and the binomial formula (A4| 



L V 



2 " 



V 2 L ^ \ n 



and obtain for (72 (^) 



n-l 



52 



L 

2n-k 



i-t) 



-k 



k=l 



1-X 



fe-r 



2n 

^ E 

k=n+l 



k = l 



L 

2n- k 



L 

+ k 



i-t) 



L 



n-k / \-x^ 



(56) 



-fc 



1 -x* 



Lv 
2i" 



— Y 



Lt 



k=l 



L 

n + k 



L 

n — k 



V. THIRD CUMULANT OF THE CURRENT 



This section is devoted to the third cumulant of the current in the stationary state of PASEP. First, we will explain 
the steps needed to derive expression (58 1 for the third cumulant, leaving the complete proof to appendix [B| We will 



then show that we recover known formulas in the case of TASEP and SSEP. Finally, we will write the expression for 
the third cumulant as a double integral in the scaling limit where the asymmetry goes to zero as the inverse of the 
square root of the size of the system. 



Exact formula 



We will now calculate the third cumulant of the current. Without loss of generality, we can suppose that L < 2n 
because of the particle-hole symmetry. It will be easier because all the sums will have their bounds between and 
2n. The steps are the same as for the two first orders: first, calculate gsit) using the recurrence relation (41 1. Then, 
insert it into (48 1 to get the eigenvalue. This time, as we do not want the next order, we will not need 63. The only 



big difference will be that these two steps use now longer expressions involving double sums (instead of simple sums 
for the second order and no sum at all for the first order). In the end of the calculation, which is detailed in appendix 
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[B] we find that the third cumulant is given by 
{L - l)Es 



p(l-x)L2 



3EE 

i>0 j>0 



6EE 

i>0 j>0 



L 

n + i 



L 

n — i 



L 

n + j 



L 

n- j 



1 + xn + x-' 

1 — X* 1 — x^ 



(58) 



L \ L \ L 

n + I J \ n + J J \ n — I — J 



n — i I \ n — J 



L 

n + i + j 



/ .9 . . .9\ 1 ~t~ X 1 ~\- X'^ 
(l^+tj+f)- 

1 — X* 1 — x-J 



3E 

i>0 



X L 

n + I J \ n — I 

L 

n 



2L\ / 3L 

l + x*y 3n{L-n)\2n) n{L-n)\3n 



1-xV 2(2L-1) / (3^-1) f L 



3 ■ 



We checked this formula numerically for systems with 2<L<12, I < n < L/2 and x — 0,0.1, ...,0.9. For all these 
cases, we computed the sequence 



f (1|(1 + 6M(7))^11) 



(59) 



for 1 < fc < 200 and e = (Xfiy(;,Tpiy i 7 being kept a formal parameter. Using convergence acceleration techniques 
(fitting our sequences with a sum of 2 exponentials using Shanks transformation), we found that, f or a ll the systems 
considered, the relative error on the third cumulant was at most 1.4 10~^, which validates formula (58 1 for the third 
cumulant. 



B. Some special cases 

1. Totally asymmetric exclusion process 

In the case of the totally asymmetric exclusion process (x = 0, the particles only hop to the right), using the 
binomial formulas ( |A3[ ), (A5) (twice), (A8l and (A9), all the sums in the expression (58) for the third cumulant can 
be calculated, giving the result 



^;3(a; = 0) 3L^n{L-n) 



2L 

2n 



4L'^n{L - n) 



3L 
3n 



P 



(i-l)(2L-l) /L 



(L-1)(3L-1) / L 



(60) 



which is the same as formula (13) in [TD]. 



2. Symmetric exclusion process 

In the case of the symmetric exclusion process (x = 1 , the particles hop to the right and to the left with the same 
rate), we see that 3 of the terms of formula (58 1 are singular. But it turns out that all these singularities cancel, as 
well as all the constants terms, giving E^, = 0. Indeed, we have the expansion near x = 1 



1 + xU + x^ 



1 — x' 1 — x^ ij (x — 1)2 ij X — 1 



Oil). 



(61) 



The coefficients of (x — 1)^^ and of (x — 1)^^ are the same, so we only have one calculation to do to cancel both the 
divergent and the constant term of Et, when x ^ 1. Using the binomial formulas (A4| and (All I, we find that both 
terms vanish, leaving us with 



E3{x^l)^0. 



(62) 
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This result was expected: all the odd cumulants vanish for the symmetric exclusion process as the system is invariant 
if we exchange the right and the left. The probability density function of the current is even. 



C. Scaling limit when a; ^ 1 as L 



The behavior of the system strongly depends on whether x = 1 or not. If a; = 1, the rates at which the particles 
hop to the left and to the right are equal and the system belongs to the universality class of the Edwards- Wilkinson 
equation. If x ^ I, the rates are not symmetric anymore and the system belongs to the universality class of the 
Kardar-Parisi-Zhang equation. The separation between these two regimes is given by the scaling 1 — a; '--^ L~^/^ with 
fixed ratio p = n/L. This scaling is a rather natural separation between weak and strong asymmetry. Indeed, a 
particle makes i? oc (1 — x)At/L revolutions through the periodic lattice during a time interval At, from the value 
of the current. The natural time interval to be considered corresponds to the time necessary for the system to reach 
its stationary state. This time scales as where the dynamical exponent z equals 3/2 for the ASEP [7j. Thus, 
1 — X ^ is the scaling corresponding to a number of rotations R ^ 1 during the time necessary for the system to 

reach its stationary state. If 1 — a; <C L~^^'^, the number of rotations during this time is <C 1 and the system behaves 
like the symmetric exclusion process. On the contrary, if 1 — x 3> L^^^^, the number of rotations is 3> 1 and the 
system behaves like the totally asymmetric exclusion process. 



In [TTj, the scaling limit 1 — x ~ l/vX of the diffusion constant was calculated from the exact expression (54 1. It 
was found (up to a factor due to the fact that they calculated the current through a bond) that 

$/i2($)p(l with /i2(*)==4/ du — -e"" (63) 

p Jq tanh(q>u) 



and 

X = (64) 

$ being held constant in this scaling limit. We will do the same here for the third cumulant in the same scaling limit. 

(66) 



From the definition (64 1 of /, we have 



l + x' 1 1 



1 tanh(^) tanh 



V-f^p(i-p)^ 

When L goes to infinity with fixed <&, / ^ 0, so that f ^ I — x and (1 — x)VL is also kept constant, which is precisely 



the wanted scaling limit. Letting L go to infinity with fixed p and the expression (58 1 of the third cumulant 
becomes, using Stirling's approximation for the binomial coefficients 



.2 , -2 



Es /p{l-p) ^v^M^ (^2 + j2)e~T7Kw) (i^ + ij+f)e' ^pC^-p) 



i>o j>o tanh ( , ) tanh ( , *>o j>o tanh , tanh -^J2= 

\^Lp(l-p) J \^Lp(l-p) J y^Lp{l-p) J \^Lp{l~p) 



j>o tanh 



^/Lp{l~p) 



We only wrote the dominant behavior of each term. The sums can now be written as Riemann integrals over the 
variables u ^ i/ yjLp{l — p) and v = j / Lp{\ — p) 

p ^'^^ ^" Jo Jo tanh($M) tanh($w) ^ ' 
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FIG. 1: Plot of the function hs{^) defined in (701, the asymmetry coefficient "I> ~ (1 — x')y L being defined in (65 1. In the limit 



where the size of the system L goes to infinity and the rates asymmetry \ — x goes to with fixed $, /i3("I>) is proportional to 
the third cumulant of the current Ez divided by "I>. The dashed line indicates the limit of /i3('I>) when $ goes to infinity, whose 
approximate value is 0.2488186. 



When L ^ CX3, we get 



P 



^^h^{^){p{l-p)f'^L''l\ 



with 



47r 



'"'*' = S75-'V„/. """" 



tanh(<i>M) tanh(<i>v) 



(69) 



(70) 



A numerical evaluation of /i3(<i>) tells us that it is a monotonic function (see fig. [T]). It grows from for <i> = to 
— Stt w 0.2488186... when ^ oo. More precisely, the behavior of h^i^) in these limits is 



A$2 

15 
167r 



3\/3 



when $ ^ 
Stt when $ — > c», 



as can be seen using 
1 



tanh <i>u tanh (f>w 
1 



1 when <f> — !■ cx) 
1 



1 /u 

tanh <i>u tanh <i>w uv^^ 3 



(- 



45v 9 



45u J 



$2 + ($4 



(71) 
(72) 

(73) 

when $ ^ (74) 



and performing the integrals in polar coordinates. 



Range of validity of the integral formula for the third cumulant 



So far, we proved the integral formula (69) for the third cumulant in the scaling limit I — x ^ l/\/L. We will now 
,ct for r 
i expref 

^^3(2:) 



see that it holds in fact for more general values of the asymmetry once we extract from it the factor 1 — x which is in 
front of the finite size expression ( 58 1 . More precisely, we will show that 



-{1-x) 



'^^^^(i „^)2^3 fQ^. i/i'- with < r < 1, 



(75) 



$ being understood as the function of x, p and L defined in equations (64 1 and (65). The case 1 — a; ^ l/^/i, which 
corresponds to taking the limit (j) — > 00, is the easiest. In this limit, all the expressions remain bounded and nonzero. 



and /i3(<&) simply converges to its (f> — > 00 limit (72) giving 
£3(1 - x» l/\/Z) 



-(1-x)(^-^)p^(1-p)^L3 when 00. 



(76) 
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FIG. 2: Plot of the skewness as of the probabihty distribution of the current as a function of the asymmetry in the partially 
asymmetric exclusion process. The rates for hoping to the right and to the left are p and xp in the scaling limit for which 
the asymmetry coefficient $ ~ (1 — x)'JTj defined in (651 is held constant when the size of the system goes to infinity. We 
represented here a^^/p/L, given in equation (791, as a function of $. The dashed curves are the asymptotics of the skewness 
when <1> ^ (as ~ <I>^) and <I> ^ oo (as ~ 



This is precisely equation (37) of [T6^ for the third cumulant when x ^ \. If we take a: = 0, we obtain the TASEP 
result [TO]. 

The limit 1 — x <C \/^/L, which corresponds to <& —> 0, is more difficult to obtain. From equation (711, we see 
that in this limit our integral expression for the third cumulant (69) is of order <^^L^/'^. We also know that the 
third cumulant is equal to when x = 1 (<I> = 0), so that is at least proportional to <&. However, it is not 
required to be of order <i>'^. Thus, if is s mall enough, subdominant terms in L proportional to may become 



formula (69 1 



larger than the ^■^i^/^ term in equation (67). In appendix [c| w e show that the correction to equation (67 1 at order 
vanish. We also compute the correction of order L^/"^ (C7) and show that it is equal, in the limit <i> ^ 0, to 



This is to be compared with the limit <f> of the leading order integral 



^3(1 - a; « l/y/L) 



pfL^ 



p 60 

If 1 — a; ^ 1 /L, these two expressions are of the same order. Thus, when 1 
is given by 

p 60 
. No other subleading correction can become larger when $ 



(77) 



a / L with fixed cr, the third cumulant 



(78) 



as each correction must be at 



which agrees with 
least of order 4>. 

We finally see that there are two interesting scaling limits for the third cumulant of the current: \ — x ^ l/\/i 
and 1 — a; ~ 1/L. Both are in fact natural separations between weak and strong asymmetry. Indeed, if one takes 
the dyna mica l exponent to be equal to 2 (symmetric exclusion process) and not 3/2 in the discussion we made before 
equation (63), we find that the natural scaling becomes 1 — a; ~ 1/L and not 1 — a; ~ l/\/i. In the scaling limit 
1 — a; ^ 1/\/L, the third cumulant is given by equation (75 1 while in the scaling limit 1 — a; ~ 1/L it is given by 
equation ( 78 ) . The first expression is valid as long as 1 — a; 1 /L 
1 - a; - 1/T 



X ^ ij L' with < r < 1. The second expression is valid for 
with r > 1/2. In the region 1/2 < r < 1, both expressions agree. 



2. Normalized third cumulant: skewness 



The skewness of the probability distribution of the current — E^/ A^^^ is given, in the scaling limit 1 — a: I/v'L, 
by (see fig. ^ 



(79) 
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In the limits 4> — > and 4> — > oo, it behaves as 





(j)3 


L 


15\/2 






L 





when $ ^ (80) 
" when (f> ^ cx), (81) 



In particular, as^/L goes to when <i> ^ 0, which corresponds to the totally asymmetric case. Indeed for TASEP 
the diffusion constant A scales as L^^"^ while the third cumulant scales as L^, so that as ~ L^/** ^ L. 

We see that, at least in the scaling limit for which the asymmetry goes to as L~^/^, the skewness of the total 
current is negative. This indicates that the left tail of the probability distribution of the total current is larger than 
its right tail, which means that the current of the system is lower than its mean value with higher probability. This 
can be understood by a simple argument [10 : in order to reduce the current of the system, one has to reduce the 
speed of only one of the particles, as particles can not overtake each other. On the contrary, in order to increase the 
current of the system, one has to speed up all the particles as the current is limited by the slowest one. It is therefore 
natural to find a negative skewness for the current. 

VI. CONCLUSION 

We have presented here an explicit method allowing the calculation of the first cumulants of the total current for 
PASEP with a finite number of particles in a finite periodic one dimensional lattice. This method, based on the 
underlying integrability of the model, relies on a perturbative resolution of the functional Bcthe Ansatz equation 
introduced in [TB] to calculate the first two cumulants. We obtained exact formulas for the three first cumulants: in 
particular, we recovered the known formulas for the first two, but with significantly less calculations, and gave an 



explicit formula (58) for the third cumulant. This formula is much more complicated than the one for the diffusion 



constant, involving five terms, two of which are double sums. However, in the limit where the asymmetry x goes to 



1 as the inverse of the square root of the size of the system, it takes a simple integral form (69 1 as does the diffusion 
constant in the same limit. 

The perturbative solution formulated here can be used to get all the cumulants of the current for the weakly 
asymmetric exclusion process, for which the asymmetry scales as the inverse of the size of the system. It gives both 
the leading and sub-leading terms (work in preparation). It would be interesting to be able to do the same when the 
asymmetry scales as the inverse of the square root of the size of the system, where nice integral formulas exist for the 
diffusion constant and the third cumulant. We also think that these methods of functional Bethe Ansatz could be 
used for other systems for which the Bethe Ansatz equations look quite similar to the ones of PASEP. In particular 
for open exclusion process, for which the Bethe equations have been found recently |25l 126] . Or for exclusion process 
with different classes of particles, for which there has already been work on the fluctuations of the current [11] . 

As the asymmetric exclusion process belongs to the universality class of the Kardar-Parisi-Zhang equation, the 
cumulants of the current for this discrete model in some scaling limit should be characteristic of the KPZ equation. 
The ratio of the square of the third cumulant with the product of the second and the fourth cumulants was studied 
numerically by Monte Carlo simulations for TASEP and for two others one dimensional growth models which also 
belong to the KPZ universality class |27j . The numerical values obtained there were then compared to the large size 
limit of the exact result for TASEP derived by Bethe Ansatz ^U\. Something similar could be done for the third 
cumulant of PASEP we derived here. 

For a random variable with gaussian probability distribution, only the first and second cumulants are nonzero. The 
skewness of a probability distribution, defined as the third cumulant, divided by the second to the power 3/2, is then 
the first sign of non gaussianity visible on the cumulants. It indicates the asymmetry of the probability distribution. 
It was recently measured for the distribution of the current in mesoscopic one dimensional systems subject to an 
electric driving field [28 . It would be interesting to find some real physical system for which the third cumulant of 
the current is given by our integral formula. 
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APPENDIX A: BINOMIAL COEFFICIENTS FORMULAS 



In the course of this paper, we used the following binomial coefficient formulas 

n-l 



i:(t)M)' = -(-irl(i) 

k=0 ^ ^ ^ ^ 



i:<t)(-i)' = -(-ir«^(^) 



E 

k=l 



L \ I ( 2L\ I ( L 



n + k I \ n - k I 2V2n/ 2 \ n 



L \ f L \ _ n{L - n) / L 



^ '"Kn + klxn — k 
k=i ^ ' 



2L \ n 



L \( L \ n(L - n) / 2L 
\n + k J \n-k J 2{2L -l)\2n 



EE 

i>0 j>0 



n + i I \n + J 



n — I — J I \ n — I 



L 
n-j 



L 

n + i + j 

1[3L 
3 V 3n 



2L 
2n 



EE(^+^-) 

i>0 j>0 



n + i J \ n + J 



L 

n-i-j 



+ 



L 

n ~ I J yn — j 



L 

n + i + j 



2 n{L - n) / L\ f 2L 

3 L \n j \ 2n 



n{L — n) 



i>0 j>0 



n + i I \n + J 



n-i-j 



+ 



L 
n — i 



n-j J \n+i+j 

_ 1 n(£ - n) / 3L 
3 3L - 1 V 3n 



2n2(L-n)2 / 



3 L2 V i 



[(„ + ;)( 



i>0 j>0 



n + j J \ n — i — j J \ n — i J \ n — j J \ n + i + j 



3 3L-1 I 3n ] 2L-1 \n )\2n ) 3 1? \ n 



EE(^+^-)K 
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'' 3Z2 I n 
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n + j 



L 
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L 

n-j 
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n + i + j 



(All) 
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L 

n + j 



L 
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L 
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L 



L 



L 



n — I J \ n — J I \ n + t + J 
j{L-2n) / L 



(A12) 



2L 



L 



n I \ n + J I \ n — J 



L 



We recall that the binomial coefficient 



is defined for an integer a as 
b{b-l)---{b-a+l) 



al 

1 for a — 
if a < 0. 



if a > 



(A13) 



In particular, if 6 is a positive integer 



if a > 6. 



Formulas (All and (A2| can be easily proved using Pascal's triangle formula 



b-1 
a- 1 



6-1 
a 



(A14) 
(once 



for the first equation, twice for the second) 



We will go on with formulas (A4|, (All 



and ( A12 1, which are the easiest. Indeed, if we call F the function of the 
re summand), one can find a "discrete primitive" G oi F with respect to 



summation index that we want to sum (t 
the summation variable, that is, if we call k the summation variable, F{k) can be written as 

F{k)^G{k + l)-G{k). 



(A15) 



Knowing this expression for F(k), it is then easy to do the summation for any range of summation. The F functions 
are respectively 



i^/c) = k 

and the G functions are 

G |^TT1 - FT2( ») = 



L 



n + k I \ n — k 



L 



L 



L 



n + I J \ n + J J \ n — I — J 



L 



n ~ I J \ n — J J \ n + I + J 



L 
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2L 



L 

n — k 
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n + I J \ n + J J \ n — I — J 



L 



{L - n + i){n + i + j) ( L 



2L 



n — I I \ n — J I \ n + I + J 



L 



(A16) 
(A17) 
(A18) 

(A19) 
(A20) 
(A21) 



which can be checked easily, and gives us the wanted formulas. These discrete primitives can be found if F is hy- 
pergeometric in the summation variable using Gosper's algorithm |29j . We used here an implementation of Gosper's 



16 



algorithm written by Peter Paule and Markus Schorn fSIT. If an hypergeometric "discrete primitive" exists, then 
Gosper's algorithm will find it. For the other single variable summations, the algorithm fails: there is no hypergeo- 
metric discrete primitive. 



We now move on to the last 2 simple sums (A3 1 and (A5l. Here, the boundaries of the sums are important: with 



other boundaries, there might not be a simple formula for the sum. We will calculate them using a suitable generating 
function. Starting from the term in t^" in the expansion in powers of i of (1 + xt)'"{l + yt)^, and keeping out of the 
sums the terms in t", we get 



E 

i>0 



i L 

n + t I \ n — t 



(A22) 



In a; = y = 1, it gives (A3 1. If we take the derivative with respect to x and to y and then set x and y to 1, we get 
( [as] ). We can now sum ( |A11[ ) and ( |A12[ ) over j and use ( [A3| and ( [A4| to prove ( |A7t and ( A10| ). Multiplying ( [All I 
by j and summing over j gives us a linear combination of ( A8 1 and ( A9 ) , namely ( AS I + j ( ^ ) . 

A t th is poi nt, o nly (A6) and either (A8l or (A9l are left. We will calculate them using a generating function, as 



for (A3 1 and ( A5 1. Starting from the term in t"^ in the expansion in powers of i of (1 + xt) (I + yt) {1 + zt) , and 



keeping out of the sums the terms in t", we get 



{l + xt)^{l + yt)^{l + zt) 
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x^ 
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y 
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L 
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L 

n- j 



L 



J 
L 

n + i + j 



(A23) 



In X 



y 



1, it gives (A6l. If we take the derivative with respect to x and to z and then set x, y and z to 1, we 



get another linear combination of (A8l and (A9l, —{A8) — {A9), proving the last two binomial formulas 



APPENDIX B: CALCULATION OF THE THIRD CUMULANT 

In this appendix, we calculate the third cumulant for systems for which L < 2n. It will be easier because all the 
sums will have their bounds between and 2n. In the end, using the particle-hole symmetry, we will notice that our 
result holds in fact for all value of L. 



1. Calculation of (73 (t) 



Noticing that for any function /, [f{xt)]^a) = x'^[f{t)]{a}i the recurrence equation (41 1 for g(t) can be written in 
the more compact form 



r-l (r-k)n , 



fc=l 6=1 



r-1 L {r-k)n kn ^c-n-^ j 

EE E E d _ jtc-a-nUbtc-l'^^ " ^"<'' " l«<c)[(l - t)^](a)[9r-k](-b)[9k](-c) 
k=l a=0 6=1 c=l ^ ' 



with 



( 1 d V , 
— i n T and 

^ if w = 



1 if M < f 

if u > V 



Introducing the notation 



(Bl) 



(B2) 



I[a,b] = [(l-i)^](a)N(-6)[5l](-„) - K-1)" ( ) N(-6), 



(B3) 
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FIG. 3: Range of summation over i and j given in equation (IBM. The + indicate that li<j+„(l — li<j — li<n) = +1, the — 
that li<j_|_„(l — li<j — li<n) = —1, and the that li<j+n(l ii<j ~ li<n) = due to a cancellation of a +1 with a —1. The 
hatched parts of the graph indicate terms that do not appear in the sums, because of li^^j, li<j+n, or because of the sums 
boundaries. The boundaries of the +, — and parts of the graph belong to these parts. The boundaries of the hatched parts 
do not belong to them. 



we get at order r — i 



2n 



L 2r, 



3 = 1 i=0 j=l 



(B4) 



which will give us the solution of equation (|8| corresponding to the stationary state provided (40) is satisfied, that is 

L 2n 



i=0 3 = 1 



(B5) 



which is easily proved. From now on, we will consider that L < 2n. Then, thanks to formula (48), we will only need 
the powers of t ranging from to in g{t) to be able to calcul ate t he third order of the eigenvalue in the case 
L < 2n. We will note this [g'i{t)\l2n- the double sum of equation (B4|, li<j-|_„(l — li<j — li<„) is either equal to 
— 1, or 1. The only indices («, j) for which li<j_|_„(l — li<j — li<n) is nonzero are (see fig. p|) 



0<i<n-l i+\<j <i + n lj<j+„(l - li<j - li<„) = -1 

n<i<2n i-n+l<j<i-l li<j+„(l - ti^j - li<„) = +1. 



(B6) 



givmg 

j — 1 i—nj—i — n-\-l i—0 

We now split the sums over j according to the expression of 92 (t) 
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We switch the sums over i and j in the double sums and then respectively make the following changes of variables in 
the sums 



Sum 1 
Sum 3 
Sum 5 
Sum 7 



i ^ n + i 

i 71 — j and i — !■ n + i 

j ^ n — j then i ^ n — i — j 



Sum 2 
Sum 4 
Sum 6 
Sum 8 



j ^ n + i 
i ^ n — i 

j n + j then i n + i + j 
j ^ n + j and i ^ n — i 



(B9) 



We obtain 
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Using the expression (57 1 for 52 (i) we can finally write, cancelling a few terms and taking the lower boundaries of all 
the sums equal to 1 
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2. Calculation of £3 



Using ([T]) and ( 48 1 , we express the third cumulant as 
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We insert the expression (Bll ) of g3{t) in the latter equation, move the sum over I inside all the other sums containing 
t and finally set I to either n — j, n + j, n — i, n + i, n — i — j or Ji + i + j according the sum using the [ ]_/. Rewriting 



1 /I 



1 — x'' 



- 1 



and 
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(B13) 
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for all A: = i, J or i + J in (Bll ) and using 



1 - V 1 - 1 - X-J 



l+x'l + x^ 
1 — x"^ 1 — X^ 



(B15) 



and the binomial formulas (A4| (two times), (A8l, (A9l and (AlOl, which allow us to calculate the sums not involving 
X, we get 
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The last two sums identically vanish thanks to the binomial formulas ( All ) and ( A12 1. This ends the proof of formula 
(58 1 for the third cumulant. 



APPENDIX C: SUBLEADING CORRECTIONS TO THE THIRD CUMULANT IN THE SCALING LIMIT 



In this appendix, we compute the two first subleading corrections to the integral formula (69) for the third cumulant 



starting from the finite size expression ( 58 1 . We will need subleading corrections to the binomial coefficient formulas 



that we used to derive the integral expression (69 1. In particular, we will use 
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and 
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f (l-3p+3p")i^ 
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1 / (1 - 2p)27;6 



£4 \72p^{l-py 
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The • • • terms contain lower powers in i that come from subdominant corrections to the Stirhng formula. They will 
not contribute to the corrections of we will calculate. In the following, we will also need an expansion of the 
difference between a Riemann sum and the corresponding limit integral. Let / be a function of 2 variables with a 
finite limit when any of the two variables goes to and decaying exponentially at infinity. We want to calculate the 
difference between the sum 



i=i j=i 



and the integral 



at order 2 in e. Writing / as 
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(C5) 



and taking care of the fact that the sums begin with i and j equal to 1, we find 

/(0,0) 



1 ^oo 1 r°° \ 
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(C6) 

Using the finite size equation (58 1, we now can calculate subleading corrections to i?3/(l — x). The term of order 
we get is the one we found previously (75). At the next order, we find that the L^/'^ term cancels out. Finally, the 
term of order is given by 
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(C7) 



This expression can be easily checked numerically: subtracting this expression plus the leading term to the finite 
size expression (58 1 for various values of <i>, n and L, we find that the remaining term is of order L ^ L^. Taking now 
the limit $ — > in expression (C7), we find that the divergent term vanishes as expected, leaving us with 



l-X 



p\l-pYL^ when$^0 



(C8) 
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